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Abstract 

We present several implementations of the Metropolis method, an adaptive 
Monte Carlo algorithm, which allow for the calculation of multi-dimensional 
integrals over arbitrary on-shell four-momentum phase space. The Metropo- 
lis technique reveals itself very suitable for the treatment of high energy pro- 
cesses in particle physics, particularly when the number of final state objects 
and of kinematic constraints on the latter gets larger. We compare the per- 
formances of the Metropolis algorithm with those of other programs widely 
used in numerical simulations. 
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1 Introduction 



Over the past few years, high energy particle physics has experienced a tremen- 
dous advance in the number of methods developed to calculate exactly scattering 
and decay amplitudes of elementary processes. Several techniques exist nowadays 
to evaluate matrix elements (MEs) for multi-particle events, both analytical and 
numerical, no matter the actual number of particles involved. 

However, in phenomenological studies, the evaluation of the MEs, in terms of the 
external particle momenta, is only a part of the whole story. In fact, in order to 
perform realistic analyses, given the finite resolution of particle detectors, a multi- 
dimensional integration over the whole or some region of the phase space available 
to the final state particles has to be performed. Indeed, in most cases, an analytical 
evaluation of the integrals become extremely complicated, if not impossible, either 
because of the complexity of the integrand and the large number of dimensions 
involved or because of the presence of cuts in the integration region]^. 

In such cases, one has to necessarily rely on numerical methods. Moreover, as 
the dimension of the integral increases, the number of evaluations of the integrand 
function needed in any generalised one-dimensional numerical approach inevitably 
grows exponentially. Therefore, the recourse to Monte Carlo (MC) methods be- 
comes mandatory: it is well known that the rate of convergence of MC algorithms 
is independent of the dimensionality of the integral . 

Naive MC algorithms are typically based on estimating the average value of the 
integrand function by sampling the latter at discrete points generated according 
to a uniform statistical (i.e., random) distribution. However, this approach turns 
out to be unsatisfactory if the integrand function strongly fluctuates over the in- 
tegration volume, as it is the case in many high energy particle processes. Thus, 
a strategy that reduces the variance of the integrand ought to be incorporated, in 
order to improve upon the tendency to converge to the correct value (see Ref. 
for details). 

Two approaches have become popular to this end, known as stratified and impor- 
tance sampling. Whereas in the latter case the algorithm concentrates where the 
integrand is largest in magnitude, in the former the function is primarily sampled 
where the contribution to the error is largest. Both techniques suffer however from 
a shortcoming: namely, the need to know beforehand the behaviour of the inte- 
grand over the all integration volume, in order to optimise the implementation. 
Unfortunately, it is exactly this knowledge that is often missing in particle physics 
phenomenology. 

^Besides, in hadron-hadron and lepton-hadron collisions, the integration over the Parton Dis- 
tribution Functions (PDFs) cannot be performed analytically, as these are accounted for by 
programs implementing numerical fits to various data sets. 
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A successful way of improving in this respect is represented by adaptive proce- 
dures. These normally involve a division of the original integration region into a 
predetermined number of 'bins' (word that we adopt here to signify any partition 
of the integration volume), with a certain number of points in each of these, the 
latter being at times redefined so to perform importance sampling automatically. 
A very much used example of such a technique is the program VEGAS 0. As a 
matter of fact, VEGAS also makes use of some stratified sampling, in order to im- 
prove the convergence in low dimensions. Because of its efficiency and accuracy, 
this algorithm has eventually come to set the standard in many particle physics 
calculations^. 

The appearance of large fluctuations in the integrand is often associated to the 
presence of 'singularities' in the MEs. In this respect, one can broadly distinguish 
between integrable singularities (e.g., resonances) inside the phase space and non- 
integrable ones (e.g., infrared emissions) at its borders. Taken separately, they 
may both be considered as 'factorisable'. In the sense that, provided a suitable 
choice (or mapping) of the integration variables is adopted, then the integrand can 
be transformed into a smoother function everywhere over the space space, with 
the MC points being generated according to a (suitable) non-uniform probability 
density. In many examples in high energy physics, however, the two kind of poles 
can occur at the same time and, possibly, there can be more than one of each type: 
particularly, as the number of final state particles, their nature and their produc- 
tion channels proliferate. Under these circumstances, the singular structure of the 
integrand becomes 'non-factorisable', in the sense that there generally exists no 
change of variables that allows even the adaptive algorithm to sample simultane- 
ously all the singularities of the integrand in an efficient manner. A multi-channel 
approach can prove useful in such cases |0. Here, the actual mapping used to 
generate a single event is chosen randomly, using a predetermined set of prob- 
abilities (or weights). A combination of the VEGAS algorithm with the adaptive 
multi-channel sampling of Ref. has recently been proposed 0. 

Another example of adaptive MC algorithm is the Metropolis technique [0 . Widely 
used in statistical physics (see Ref. |§] for a review), it has nonetheless seen very 
little applications to particle physics. These have been mainly confined to the case 
of Lattice Gauge Theories (LGTs) . It is our intention here to demonstrate its 
potential in the context of a Quantum Field Theory (QFT) of the continuum: in 
evaluating total and differential cross sections in multi-particle scattering processes 
at high energies. In particular, we will show how a Metropolis algorithm easily lends 
itself to several manipulations, that make it a versatile instrument in performing 
such calculations, thus overcoming most of the problems that we have described, 
often matching in speed, efficiency and accuracy VEGAS itself and outclassing many 
others of the most commonly used algorithms. 

^For parallel versions of VEGAS, see Ref. 
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The plan of the paper is as follows. In Sect. |^ we describe the fundamentals 
of the Metropolis technique and propose several integration methods. Sect. ^ 
presents a few numerical examples and comparisons with other algorithms in the 
context of some benchmark processes in high energy particle physics. In Sect, ^we 
summarise our results. Finally, we will defer to the Appendix the analytical proof 
of a condition required to the algorithm for its convergence in four-momentum 
phase space. 

2 The Metropolis algorithm for a QFT 

The Metropolis method is a somewhat adaptive MC integration algorithm which 
is widely used in numerical statistical mechanics and LGTs. We will here demon- 
strate how it can be modified to evaluate cross sections and other observables from 
perturbative ME calculations of scattering and decay processes in a QFT. The 
integration is performed over the four-momentum phase space of the final state 
particles. The phase space can be of arbitrary form, with the only requirement 
that the particles are on their mass-shell, i.e., = mf = constant, where pi and 
rrii represent the four-momentum and mass of the i-th particle. This restriction is 
enforced by construction, in order to give a correct description of the phase space 
associated to the final state particles. However, intermediate particles, that can 
appear in a process, are allowed to be off-shell. 

2.1 Description of the algorithm 

In general, the problem is defined by a phase space and a weight-function 
for every point x in phase space. A random walk in the phase space is performed 
by starting at an arbitrary initial point xq, and generating new points Xi by using 
the weight w. In the Metropolis algorithm, the sequence {xi} is generated in a 
way which ensures that the points will reach the correct distribution when the 
number of steps (hereafter, A^) is large. In a QFT, the phase space is given by 
the on-shell four-momenta of the outgoing (and possibly incoming) particles, the 
weight-function coincides with the ME, and the cross section will correspond to 
the partition function (defined below) of the problem. 

In it is shown that in order to reach the correct probability distribution, P{x) oc 
w{x)dx, for the sequence, one has to fulfill two conditions in the stepping procedure 

Xi ^ X^-j-i. 

•^In statistical mechanics, w is normally given by an energy-function, a temperature and the 
Boltzmann distribution. 
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(a) All points in phase space must be reachable with a finite number of steps. 

(b) The condition of 'detailed balance' must be fulfilled: 

P{Xi)P{Xi Xi+i) = P{Xi+i)P{Xi+i Xi). 

The normal procedure to satisfy detailed balance (b) is to randomly choose a point 
Xi+i with even distribution inside a region fij. This new point is then accepted 
with probability P = mm{l,w{xi+i)/w{xi)). If it is not accepted, the previous 
point will be put in the sequence once more: i.e., Xj+i = Xi. Condition (a) is then 
satisfied if the overlap of {fli} can cover the whole of the phase space. Choosing 
the fij's to actually be smaller than the phase space will make the Metropohs 
algorithm adaptive since, at each step, the suggested point Xi+i is likely to be in 
the region of large weights. 

For any observable, 0{x), the mean can then be calculated as follows: 



<o>. 



J 0{x)w{x)dx 
J w{x)dx 



— lim 



N 



(1) 



The partition function, Z, can be calculated by evaluating the average of l/w and 
using the relation: 

Z = [ wix)dx = ^^"^ . (2) 

However, a straightforward use of the Metropolis algorithm will in this case be 
inefficient since 1 / tv has large contributions when w is small. The regions of small 
w are not visited so often, and one would have large statistical fiuctuations. For a 
w that has large variations, it is then better to weight the random walk with w^/'^ 
instead, and to evaluate Z by: 



Z = ( w(x)dx — . ,n ^ ^ f dx. 

J ^ ' <W-l/2>,„i/2 i 



(3) 



In this way, the magnitude of the variations in the weight-function are effectively 
halved. The total volume of phase space, / rfx, has to be estimated separately. 

In a QFT one is often interested in evaluating the total cross section other than 
observables of the four-momenta. For example, for processes with two particles in 
the initial and n in the final state (scattering reactions), the cross section can be 
written as 



a — 



2pL 



/n 
\M\\Va.V^AVr'i)Vi 



) S{ptot-Y.Pi) ' (4) 



where ptot is the total four-momentum, Pa and Pb are the incoming momenta, {pi} 
are the n final state momenta, with rrii their respective masses, and |A^p the ME 
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associated to the reaction considered. By choosing the weight-function and the 
phase space to be 



(27r 



\4-3n 



W = \M\ (Pa,Pfe, {Pi}), (5) 

'^Ptot 

n 

dx = Yi [d^Pi^^ {pI -^1)]^ [ptot -^Pi) , (6) 

i=l 

respectively, we can use the Metropohs algorithm to perform the integration. 

In order to do so, one has to guarantee that the stepping procedure in the Lorentz- 
invariant phase space can be done in a way that satisfies the two conditions 
(a) and (b) above. We now describe how this can be achieved: first comes an 
illustration (points (i) to (iv)) of how the suggested point in phase space, Xj+i, can 
be generated; afterwards, we show that this procedure complies with (a) and (b). 

The proposed point in phase space, Xj+i, can be generated in the following way. 

(i) Choose two of the final state particles randomly. 

(ii) Boost them into their centre-of-mass (CM) frame. 

(iii) In that frame, rotate them randomly with an even distribution in 47r. 

(iv) Boost them back into the original frame. 

In order to demonstrate that this procedure will satisfy conditions (a) and (b), we 
will make use of the Lorentz-invariance of dx in eq. (P). To make sure that (b) 
is fulfilled, we have to show that the suggested point is generated with an even 
distribution inside the region fi, which is reachable in one step. But it is clear 
that this is true in the CM frame of the two particles chosen in (i). Therefore, and 
because of Lorentz-invariance, it is true in any frame. 

To demonstrate that also condition (a) is fulfilled, we have to show that from any 
point in phase space, {pi}-, one can reach any other point, {p-}, by a finite number 
of steps. This can actually be achieved in, at most, n — 1 steps. Let Apj be the 
shift of the i-th four-momentum, p'i= Pi + Apj. Choose an arbitrary particle, say 
Pn- For each of the n — 1 other momenta, Apj can be transferred to p„ in one step. 
This is clear since pi + pn is conserved. In the CM frame of pi and p„, the whole 
surface of conserved pi + pn and the other momenta fixed can be reached in one 
step. So is in any frame: any Apj which conserves Pi -f- p„ can be transferred in 
one step. By doing this for all the n — 1 other particles, we have reached the point 

In some cases, one might want to introduce cuts in the integration volume^. That 
is, instead of integrating over the whole of phase space, PS, one might want to 

Alternatively, the ME might be identically zero inside a finite region of it. 
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integrate inside a reduced region, P5'red- In that case, one has to make sure that 
condition (a) is fulfilled also over PSj-ed- In the Appendix we show that this is 
the case if PS^ed is a connected region. If this is not true, then PSred has to be 
separated into parts, PSred,i, each of which is connected, and the integration be 
done for each of these separately. 

Finally, in some cases, it can be desirable to integrate a ME (or just a term of it, in 
which case we use the notation T) which is negative over some parts of the phase 
space. If so, T can not directly be used as the weight-function w. Instead, one can 
use its absolute value, and keep track of how often T is negative. With w = |T|, 
the partition function is in this case given by: 



In summary, the average of an observable O, given by a function 0{{pi}), can be 
calculated as follows. 

• Choose a suitable starting point {pio}. 

• Generate a new point {pn} as described above. 

• Accept this point with the probability 



with w given by the ME (eq. (||)). 

• If not accepted, use the previous point: {pn} = {pm}- 

• Store the value of O: Osum = C'sum + C^iiPn})- 

• Repeat the four last steps until the desired level of convergence is reached 
{N steps). 

• Finally, take the average: <(9>~ 

2.2 Time consumption and error estimates 

In most of the (sub)processes that we have integrated with the Metropolis al- 
gorithm, the, by far, most time-consuming tasks have been the calls to the ME 
function. Our strategy to reduce the CPU-time has then been to, as far as possible, 
reduce the error for a fixed number of calls to this function. 




(7) 
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The algorithm requires many boosts and rotations of four-vectors. For 1 mil- 
hon steps through phase space, e.g., this takes around 30 seconds on a 350 MHz 
Pentium-II processor, on a LINUX platform. Calling the ME the same number of 
times on the same machine is often a matter of hours. Thus, the stepping proce- 
dure has not been optimised. Instead, we have tried different ways of refining the 
interface between the algorithm itself and its use of the ME. 

A complication in evaluating the statistical error in a Metropolis-based integration 
is that the generated points in phase space are correlated. The correlation length 
effectively reduces the number of statistically independent data points. We have 
handled this flaw by collecting the data into a number of 'blocks'. The average in 
each block is then taken and used as a new, statistically independent data point. 
The size of the blocks must be made larger than the correlation length, while the 
number of blocks must be large enough, so that the variance evaluated from the 
averages can be trusted. The normal procedure to check that the correlation length 
has been reached is to divide each blocks into smaller parts, and check that single 
block-parts inside one block are not more correlated than the block-parts from 
different blocks, see Ref. [|T^. 

2.3 Some methods of refinement 

In this Subsection, we will describe some methods of refinement which are devoted 
to reduce the error in the integration for a certain number of calls to the ME. These 
methods will in general increase the CPU-time of the algorithm significantly, but, 
as mentioned above, the overall CPU-consumption will increase marginally, while 
improving the efficiency of the algorithm in various respects. 

2.3.1 The parallel integration method 

In realistic problems, the ME can present large variations, often due to resonances 
arising from intermediate particles. In simple terms, large variations make it hard 
to move around in phase space. The correlation length in the Metropolis steps will 
be large and many calls to the ME[| will be required in order to get independent 
data points. We will here describe a method, the 'parallel integration method', 
which enables one to move large distances in phase space, without using the exact 
(time consuming) ME function. 

First, one can introduce an approximate ME, |A1p(x). This should be a simple 
function which is quick to call, and which is a good approximation in the regions 
where has large contributions. It could, e.g., be a product of just the 

^Hereafter, denoted by the short-hand notation |A^p(2;) = \A4\'^{pa,Pb, {pi}), where x repre- 
sents the phase space point {pt} of eq. (^. 



7 



resonant propagators in |A1p(a;). Then the configuration space is enlarged with 
an extra, discrete parameter, s, which can assume the values or 1. A weight- 
function w is constructed in the enlarged configuration space, such that it returns 
\A4\'^{x) if s = 1 and |A^p(a;) if s = 0. Finally, the Metropohs-step is modified so 
that, occasionally, instead of suggesting a new set of four-momenta, a new value of 
s is suggested. Data points are to be taken only when s = 1. The distributions of 
points in the part with s = 1 (of our enlarged configuration-space) will be exactly 
as in the original phase space but with only |A^p(x) as the weight function. We 
know this since, in general, the Metropolis algorithm ensures that the generated 
points are distributed according to the weight-function, and this is of course true 
also for individual parts of the configuration space. 

The gain of this method is that, during the steps in phase space when s = 0, one 
can reach distant points, without calling the slow, exact |A1p(x). At will, one can 
choose to make it more probable that s = 0, so that |7Mp(a;) is used as seldom as 
possible. For example, this can be done by introducing two integers, A^^o and A^i. If 
s = i, we attempt to switch the value of s after Aj steps. We can then choose, e.g., 
A"o = 100 and A^i = 1. The optimal value of Nq/Ni will be given by how faster the 
approximate ME function actually is in comparison to the exact one. 

Let us summarise the parallel integration method. 



Find a crude and quick approximation, to the original ME, \ Ai 

Enlarge the phase space by introducing the discrete parameter s = {0, 1}. 
Construct the weight function 



x] 



w{x, s) 



\M\\x) if s = 1, 
\MWx) if s = 0. 



Choose an arbitrary starting point Xq in phase space and set s = OQ. 

Perform Aq Metropolis steps with the weight w{x, 0) = and without 

collecting the value of the observable O. 

Switch the value of s with probability 

P = mm 1, , , = mm 1, ' ^ ^ ' 



'w{x,0)J \ ' \M\^{x)J 

If the switch is accepted, store the value of the observable. 



'Also arbitrary, but s = is recommended. 
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• The last three points are repeated in the following way until desired conver- 
gence is reached: 

— li s = i, perform Ni Metropolis steps and, only if s = 1, store the value 
of O after each step. 

— Switch the value of s with the probability 

„ / w(x,s')\ 

P = min 1,^^ , 

where s' = 1 if s = and s' = if s = 1. 

— Store the value of O, if the new s = 1. 

• Take the average of the observable. 



With this method, the total cross-section is calculable in a more effective way, 
than described before. This is the case if the cross-section of the approximate ME 
is already known with high accuracy. By counting the number of times, rj, that 
s = i after a switch, one gets the relative magnitudes of the cross-sections from the 
fraction of the r^'s. If the exact cross-section is denoted by a, and the approximate 
one by a, one has that 

a = a-. (8) 



2.3.2 The variable energy method 



In some cases, one is interested in doing the integration with a dynamically variable 
total energy, i?tot 7^ constant. This might be the case, e.g., if one is interested in 
the cross section as a function of i^tot in a certain range, or if the two incoming 
particles can have varying energy, e.g., depending on some given PDF. To this 
end, we describe here the 'variable energy method', where the configuration space 
is extended by adding i^tot as a dynamical parameter. Allowing for a variable 
energy can actually make the integration more effective, for the following reason: 
it is often the case that the variations of the ME are larger for larger energies. 
The phase space walk is thus easier at low energies. If one is interested in some 
observables to be calculated at large -Etot, it can be easier to reach farther in phase 
space by, say, taking a round tour into the lower energy regions. This idea has been 
explored in calculations in statistical mechanics, where often is the temperature to 



be used as a dynamical variable 12 



One proceeds as follows. Let a point in the extended configuration space be denoted 
hj X = {{pi},Etot), where £^tot has been written explicitly, though it is actually 
computed from the four-momenta. In order to be able to dynamically take steps 
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to other energies, we introduce a one-to-one mapping {pi} = K {{pi} , Et^t) , which 
gives a new point, {{pi}, -Etot), given the previous one, {{pi}, -Etot), and the new en- 
ergy as well, Etot- We also need a phase space weight, p[{{pi}, -Etot) ^ {{Pi}, -^tot)], 
which is evaluated from phase space densities at the different energies (as described 
below). In the Metropolis path, one can alternatively attempt to change the energy 
or just update the configuration. A step which can change i^tot is then taken in 
the following way. 

• Suggest a new energy, Etot, chosen with even distribution in a certain region. 

• Find the corresponding point in phase space, {pi} = K {{pi} , E^ot) ■ 

• Accept this with probability 

\Mmp.})\ 



P = min l,p[{{p,},Etot) ^ (fe},^tot)] 



How the mapping K can be defined and how the corresponding phase space weight 
p is calculated is the next step. In order to describe how this can be achieved, we 
have to digress briefly, by giving a general description of the Lorentz phase space, 
defined in eq. (|^) (for an overview, see, e.g., Ref. pTSl). 



Let Vn{{'mi}, s) denote the total volume of the n-particle phase space with masses 
{rrii} at the squared invariant energy s. For = 2, we have 

V2{ml, ml, s) = —\J\{ml, ml, s), (9) 

with X{a, b, c) = a"^ + b"^ + — 2ab — 2bc — 2ca. In the CM frame, the magnitude 
of the two outgoing momenta is given by 

P{ml, ml, s) = ^^X{ml,mls). (10) 

The volume for n particles can then be calculated recursively: 

Vn{{mi}n,s)= ds'V2{s',ml,s)Vn-i{{mi}n-i,s'), (11) 



where the integration range is restricted between sq = (j2i ^ t^i) and si = 

{^/s — mn) . The mapping K, from the energy E^ox. to -Etot = E^ot + ^E, can 
be defined by changing the momentum of the n-th particle and boosting the oth- 
ers, so that the same CM frame is kept. In case this is actually possible, the 
back-to-back momentum of the n-th particle and the n — 1 particle system, in the 
CM frame, is given by P{s' ,m'^, E^^^), with the momentum-function P as in eq. 

(0)- 
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When the change in energy can be done in this way (i.e., provided A > 0), p is 
given by the ratio of the 2-particle volume for the two energies. This is the case 
since the value of s' is preserved and the n — 1 particle volume in the integrand in 
eq. ([TI|) is not affected. Thus, the p expression corresponding to our mapping K 
is given by 

where the ^-function returns a zero when A < (corresponding to the case in which 
s' is out of the integration range of eq. (0)), and 1 otherwise. We stress again 
that the choice of the mapping K is not unique. We have illustrated here a simple 
method, however, it is possible that other, more sophisticated mappings can be 
chosen to make the whole procedure more efficient. 

2.3.3 Variable energy and parallel integration 

The two methods described above can be used simultaneously, as follows. Choose 
an energy interval IE, where the integration is to be performed, and a fixed energy, 
Ep, to be used only for the approximate ME function |A1p(a;). Also define a smaller 
interval IE^/ C IE, which will be the 'window' where we switch between the exact 
and approximate MEs. The need for such a region is induced by the fact that the 
cross section often varies a lot with energy. The window should be put around 
the region in energy where the cross section is expected to be large. It is also 
recommendable to have the fixed energy, Ep, inside this window. The algorithm 
then goes as for the parallel integration method except that, when the exact ME 
is to be used, one allows also for energy updates. We here describe how the cross 
section as a function of energy is evaluated, assuming that the cross section for the 
approximate ME, a{Ep), is known at one energy, Ep. 

The stepping procedure is done as sketched below. 

• Choose the energy ranges IE D lEw, and the point Ep G lEw- IE is the 
range we are interested in, lEw should be chosen in the region of large cross 
sections and Ep is the point where a is known. Let, as before, i^tot denote 
the latest accepted energy used for the exact ME. The starting value could, 
e.g., be Etot = Ep. 

• Introduce the discrete parameter, s = {0, 1}, and the two update periods, 
Nq and Ni, as before. 

• Start with s = at the energy Ep and perform Nq ordinary Metropolis steps 
with the approximate ME. 



11 



• Attempt a switch of the value of s by proposing an energy step Ep — > E^ot, 
as described in Subsect. p.3.2| , with the new weight given by the exact ME. 

• Whenever s = 1 is accepted, an energy-varying Metropohs path is taken in 
the region IE. Whenever A^^i points inside the window IE\y have been chosen, 
a new switch of s is attempted to the point Ep and with the approximate 
ME as the new weight. 

In order to get an estimate for the total cross section for the exact ME, cr{E), one 
needs to store the distribution of energies generated when s = 1. This distribution 
is then normalised by means of a{Ep), together with an evaluation of how often 
s = 1. That is, proceed as follows. 

• Whenever s = 1 and an energy step inside IE is taken, put the newQ energy 
into a distribution, f{E). 

• Evaluate the number of times, r^, that the value of s = after a switch 
attempt. 

• Do this until desired convergence of the distribution, f{E), and the numbers, 
Tj, is reached. 

The cross section then becomes, cr{E) = Nf{E) with the normalisation factor 
= a{Ep)ri/ fro, where / is the average of f{E) inside the window lE^. 



3 Examples and performances 

In this Section we report about some numerical results obtained by using the 
Metropolis algorithm and compare them to the outputs of other integration pro- 
grams widely used in particle physics calculations. In Subsect. |0| we study multi- 
photon production in electron-positron annihilations, as an example of the ability 
of the various algorithms to deal with the problem of singularities associated with 
an increasing number of infrared emissions of different topology and over an increas- 



ingly reduced phase space. The benchmark process considered in Subsect. 2]2]is the 
production of 6-quark and VT^-boson pairs at future lepton colliders, in order to 
test the performances of the algorithms in presence of massive particles and diver- 
gent/resonant invariant mass poles, the latter further occurring in a multi-channel 



environment where use of negative weights is made. Finally, in Subsect. we 
will focus our attention to the case of some radiative top-antitop events as pro- 
duced in hadron-hadron collisions, leading to eight-parton final states, the latter 



^ This could be the same energy as before, in case the suggested energy is not accepted. 
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integrated over a reduced phase space: here, we account for the behaviour of the 
integrator as induced by a non-fixed partonic energy, a large number of particles 
in the final state and the enforcement of experimental cuts through 6'-functions in 
the integrand. 

We have carried out our tests on several machines, without finding any significant 
dependence of the relative algorithm performances upon a particular choice among 
them. For reference, we list here the platforms used: a-DECstations 3000 Model 
300 and 400 and a-Servers lOOOA 5/300 and 1000 4/200 (running both VMS and 
OSF operating systems), UNIX Hewlett-Packard Workstations type A 9000/712 
and A 9000/782 and the LINUX system already mentioned. 



3.1 Multiple infrared radiation 



The reference paper for this Subsection is |Tj] . There, it was studied the tree- level 
process e^e~ — > n 7 (with massless electrons/positrons), where 7 refers to a photon 
with n taken up to 7, and where the relevant MEs were presented analjd;ically. In 
our forthcoming tests, we have made use of the expression for the latter given in 



eq. (17) of ||T4|. Notice that such a formula is exact for n < A only, whereas it is an 
approximation for n > 4. However, being much faster than the exact ME, eq. (16) 
of that paper, while retaining its main dynamic features, it is much more useful to 
our purposes. 

The challenge here is to integrate a ME that diverges when a photon becomes either 
soft or collinear (with one of the initial state fermions). To render the results finite, 
while still allowing for the infrared behaviour at the edges of the (reduced) phase 
space, we impose the following cuts (same as in Ref. |T^ ): 

E^>5GeV, cos^^<0.9, (12) 

on the energy and (cosine of the) polar angle of each final state photon. For 
comparison purposes, we also adopt the same CM energy used in [^, that is, 
y/s = 100 GeV. However, we will increase here the number n of photons produced 
to 9. 

The algorithms used for this example, other than the Metropolis one, are VEGAS, 
RAMBO and the NAGLIB subroutines DOIEAF and DOIGCF. Of VEGAS, we have 



already given a description. RAMBO [15] is not exactly an integrator, though it 



can be used in some instances in such a fashion. Rather, it is a multi-particle 
phase space generator, as - given the total CM energy and the number of particles 
required with their masses - it produces a set of four-momenta and the phase space 
weight associated with that configuration. However, when the integrand function 
does not present sharp peaks (as it is the case here), it can be used to estimate the 
total cross section and its standard deviation simply as an arithmetic average and 
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through a standard quadrature formula, respectively. Indeed, this is the algorithm 
that was successfully employed in Ref. [Q. The subroutines DOIEAF and DOIGCF 
are part of the NAGLIB library |^. They both are mult i- dimensional adaptive 
quadrature integrators, the first over an hyper-rectangle whereas the second on a 



general product region. DOIEAF makes use of the algorithm described in Ref. ||17 



whereas for DOIGCF one should refer to |18]. Their usage and characteristics are 



well introduced in Ref. |]T6[, so we do not dwell any longer on them. The Metropolis 
implementation used here is the one described in Subsect. |2]l|, making use of the 
more effective formula in eq. (|^). The reduced phase space volume (after cuts), 
i.e., / dx in (^, was calculated numerically but with insignificantly small errors. 
Alternatively, one can assume it to be given as an exact input, if the phase space 
integration can be performed analytically. 

In this test, we have proceeded as follows. For a start, we have fixed the number 
of MC points generated to carry out the integration, A'mc; to be approximately 
10^ (including those eventually rejected because of the cuts) for all algorithms 
and irrespectively of the number of photons generated. One may consider this 
condition as a prototype of what inevitably occurs in numerical computations, 
when one can only dispose of a finite amount of CPU (here, corresponding to the 
time needed to evaluate 10^ times the integrand function, as it would actually 
happen if all momentum configurations generated were accepted)^. Notice that 
we have introduced the cuts through theta functions in energy and angles, while 
maintaining as upper and lower limits of the integration variables chosen those 
needed to cover the all of the original (massless) n-particle phase space|. Under 
these conditions, what one would expect from an optimal algorithm is both the 
tendency of promptly adapting itself to an increasingly reduced integration volume 
- as n gets larger - (thus minimising the loss of MC points through the cuts) 
and the ability to efficiently dispose of those points that survive the kinematic 
constraints (by minimising the error of the integration). A measure of the former 
is the ratio between accepted and generated MC points, whereas for the latter of 
the percentage spread of the errors about the central values obtained. We will see 
that the estimates will roughly be consistent with each other for all algorithms up 
to n = 8, while some of the programs will fail to converge for n = 9. Then, for 
those programs that manifest the ability to converge to the correct value of the 
cross sections with good accuracy up to n = 9, we have increased Amc by ten and 
hundred times, this way further checking the rate of convergence of the algorithms 



®In fact, here the expression of the rt-photon ME is considerably simple that the time spent 
in evaluating it has little impact on the total one employed by the all integration procedure, 
irrespectively of the value of n. 

^For consistency, we have customised the choice of the latter to be the same invariant masses 
and angles for VEGAS, DOIEAF and DOIGCF. As for RAMBO, there is no need to provide the integration 
variables and corresponding Jacobian factors, as this is done automatically by the program: see 
Ref. pa for specific details. 
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when, ideally, an infinite CPU (i.e., number of MC points) is made available. For 
sake of illustration, we have reproduced in this case the e~^e~ -^77 cross section. 
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Figure 1: Acceptance rates. A, defined as the ratio of accepted over generated MC 
points (times 100), for the VEGAS algorithm in evaluating the total cross section 
for e^e~ — n 7, with n = 1, ...9, for five possible choices of ITMX and NCALL, 
such that their product is approximately equal to A^mc = 10^. 



Before proceeding to compare the algorithms, one should recall that in VEGAS 
there exist two parameters, NCALL and ITMX, that determine the actual number 
of points generated. The first is the (approximate) number of integrand evalu- 
ations per iteration whereas the second is the maximum number of the latterQ: 
to change one or the other affects the overall performance of the algorithm in 
various respects [§. Thus, as a preliminary exercise, we have run VEGAS vary- 
ing these two parameters: e.g., by taking ITMX = 1,2,4,5, 10 and, consequently, 
NCALL = A^Mc/ITMX, being Nmc = 10^ The outcome is presented in Fig. |, 
where we show the acceptance rate, A. There, one can first notice the indifference 
of VEGAS to the choice of NCALL and ITMX for n < 4, that is, if the dimension of 

^°As we set the required accuracy to be negative, i.e., ACC < 0, all of them are performed. 
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the integralQ, NDIM, is below 7. For n > 5, or equivalently NDIM > 10, if one in- 
creases ITMX, the adaptabihty worsen for high dimensionahties while it improves 
for low ones. As a compromise between the two tendencies, we will adopt for the 
reminder of this test the choice ITMX = 4 and NCALL = A^mc/ITMX for any n. 
In the case a VEGAS iteration fails to find points above the cuts (as it can happen 
for very large n), thus yielding a zero, its contribution to the total estimate of the 
cross section is discarded altogether. 



(7(e^ 



n 7) (nb) at ^ = 100 GeV 



n 


Metropolis 


VEGAS 


RAMBO 




2 


2.66782 ±0.0029 


2.664973 ± 0.000050 


2.66114 ±0.0023 


xlO-2 


3 


6.34056 ± 0.029 


6.26009 ±0.017 


6.26585 ±0.014 


xlO-^ 


4 


6.63656 ± 0.036 


6.63202 ±0.050 


6.67522 ±0.025 




5 


4.0304 ±0.034 


3.94629 ±0.054 


3.99786 ±0.023 


xlO"*^ 


6 


1.49256 ±0.013 


1.48609 ±0.048 


1.47493 ±0.010 


X 10-10 


7 


3.72003 ± 0.037 


3.73682 ±0.26 


3.63071 ± 0.040 


X 10-13 


8 


5.92104 ±0.072 


5.51518 ±0.51 


5.83508 ±0.077 


xlO-i^ 


9 


6.43459 ± 0.079 


5.62617 ± 1.3 


6.57829 ±0.21 


X 10-19 



n 


DOIEAF 


DOIGCF 




2 


2.66510 ±0.0019 


2.66493 ± 0.000010 


xlO-2 


3 


6.19822 ±0.29 


6.15469 ±0.031 


xlO-^ 


4 


6.47451 ±4.1 


6.60557 ±0.13 


xlO-6 


5 


5.13410 ±4.8 


4.11628 ±0.45 


xlO-s 


6 


11.2745 ±42. 


1.44837 ±0.21 


X 10-10 


7 


0.776428 ± 33. 


2.87767 ±0.64 


X 10-13 


8 


3.56037 ±385. 


4.93633 ±2.2 


xlO-i*^ 


9 


59.0009 ±866. 


cannot compute 


X 10-19 


E^>5GeV cos^^<0.9 



Table 1: The cross sections and relative errors for e^e — >■ n 7 at ^/s = 100 GeV, 
with n = 1, ...9, as obtained from the five algorithms documented in the text. For 
VEGAS, the setup ITMX = 4 and NCALL = 250000 was adopted. Several values 
of MINCLS (see Ref. ||16|) were used for DOIEAF, but no significant improvement 
was found compared to the data reported. 



"'^"'^Note that we have mapped the n-photon phase space in such a way the the azimuthal angle 
around the electron/positron beam direction is one of the integration variables. Being the cross 
section independent of this variable, the phase space integral can be reduced by one dimension 
and simply multiplied by 2tt. Of course, four of the initial dimensions of the phase space 
integral are removed by the (5-functions associated to the four-momentum conservation between 
the initial and final states. 
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In Tab. |T| we present, as function of n, the central values and the associated er- 
rors produced by the five integrators considered in evaluating the cross sections for 
e^e~ — *• n 7 (compare to Tab. 1 of Ref. [|I^] for n < 7). (DOIGCF cannot compute 
integrals with more then 20 dimension, so that the cross section corresponding to 
n = 9 does not appear in the table.) A first obvious result (apart from the short- 
comings of DOIEAF as n increases) is that whereas VEGAS performs undoubtedly 
better than Metropolis for small n, say, below 4, if n > 4, Metropolis yields a 
much more accurate answer. Even RAMBO, a non-adaptive algorithm, excels over 
VEGAS for large photon multiplicities, though, for n = 9, not as well as Metropolis. 
The error from DOIGCF is significantly larger than that for the other algorithms (ex- 
cept DOIEAF) for n > 4, whereas for smaller values it almost achieves the accuracy 
of VEGAS. 



10-^ 



10-= 



10^ 



10^ 



10 



10 



W 10 



10 



-1 



-2 



-3 



-4 



n I I I I I 1 1 r 

O 

□ 



10'^ 



10^ 



e e ^7 photons 



n I I I I I 



□ 



10' 



□ 



10' 



<> Metropolis 
□ VEGAS 
+ RAMBO 



_l I I I I I I I I I I I I I I I L 



10^ 



10 



10 



Figure 2: Acceptance rates, A (above), defined as the ratio of accepted over gen- 
erated MC points (times 100), and relative error, E (below), defined as the ratio 
of the standard deviation over the average value, for the Metropolis, VEGAS and 
RAMBO algorithms in evaluating the total cross section for e+e^ ^ 7 7, as a function 
of A^MC- For VEGAS, the setup ITMX = 4 and NCALL = A^mc/ITMX was used. 
(Note the overlapping errors for Metropolis and RAMBO.) 
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Therefore, for high-dimensionality phase spaces, Metropohs, and RAMBO as well, 
seem to be much more accurate than VEGAS. However, one might well wonder what 
is the actual number of points used in the evaluation of the integral, as the accuracy 
of the latter strongly depends upon it. It turns out that the acceptance rate of 
VEGAS (recall Fig. I) is very poor compared to that of Metropolis and RAMBO, which 
is about 61(59)[56]% and 19(10)[4]% for n = 7(8)[9], respectively. (In particular, 
the large difference between the acceptance rates of Metropolis and RAMBO for the 
case n = 9 should explain the much smaller error for the former.) Thus, it is 
not surprising to see a bigger error in the former. Indeed, if VEGAS itself is run 
in non-adaptive mode (i.e., ITMX = 1) its acceptance significantly improves (see 
again Fig. P and the error consequently diminishes (typically halved), still being 
larger than in Metropolis and RAMBO, though. 

However, in order to show that the higher accuracy for Metropolis, as compared 
to VEGAS, for a large number of dimensions is not an artifact due the specific 
value adopted for A^mc, we plot in Fig. ^ both the acceptance and the size of the 
relative error for the two algorithms in calculating, e.g., a{e^e~ — >■ 7 7), with 
Nmc = 10^10^ and 10^ (as usual, ITMX = 4 in VEGAS). Similarly, we proceed 
for RAMBO. One can see that, no matter how many MC points one can dispose of 
in VEGAS, both the adaptability and accuracy in Metropolis remain significantly 
better. The gain for VEGAS is in the end only appreciable against RAMBO: Metropolis 
still stands out as the best choice for large n values, whatever Nmc is actually used 
for the integration. This statement remains true for any choice ITMX = 1, 2, 4, 5, 10 
adopted in VEGAS. 



3.2 Mass singularities, mult i- channels and negative weights 



The physics concerned with our discussion below can be found in Refs. [|^,^ . The 



process calculated is e~^e~ — >• bbW~^W~, with massive quarks (and gauge bosons, 
of course), proceeding at tree-level through the 61 Feynman diagrams depicted in 
Fig. 1 of Refs. ||19|,^ (again, we assume rrie = 0). These include several graph 
subsets (eight of these were isolated in Ref. pO|), each having a peculiar (non- 
) resonant structure, so that they can be regarded as separate production modes of 
an actual multi-channel process. Furthermore, since interference terms exist in the 
full ME among the various channels, some of the latter can give rise to negative 
contributions in the integration procedure (the T weights of Metropolis discussed 
previously) . 

For sake of illustration and comparison among the algorithms, rather than gen- 
erating the total e^e~ bbW^W~ cross section in a unique run using a-phori 
weights to choose among the various channels (i.e., a la [^), we instead perform 
a separate integration over each of the latter, as they present different challenges 
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to the algorithms, see Refs. |T^,|2D|. In general, our approach can be viewed as 
a preliminary by-hand optimisation of the weights eventually used in a full ME 
multi-channel run, rather than the automatic one discussed in Ref. |Q. Numeri- 
cal values used here for the various parameters needed for the calculation are as 
follows (the reader should not mind their obsolescence, as they are used for illus- 
tration purposes): = 300 GeV, rut = 145 GeV, Tt = 0.78 GeV, rub = 5 GeV, 
Mh = 120 GeV, Th = 6.9 MeV, Mz = 91.1 GeV, Tz = 2.5 GeV, Mw = 80 GeV 
and Tw = 2.2 GeV. Note that we do not impose any cuts on the phase space (so 
that all A^MC generated points are actually used in the ME evaluations) and we 
neglect initial state radiation (ISR), thus identifying the partonic energy with that 
of the collider. 



a(e+e- ^ bbW+W-) (fb) at = 300 GeV 


channel 


Metropolis 


VEGAS 


RAMBO 






598.211 ±44. 


582.541 ±0.26 


517.794 ±56. 






2.94538 ±0.12 


2.78886 ±0.019 


2.66227 ±0.64 




% 


2.85556 ±0.14 


2.79228 ±0.019 


3.82429 ±0.61 




% 


1.96508 ±0.041 


1.85213 ±0.0013 


1.82016 ±0.082 


xlO-2 


% 


4.99161 ±0.079 


5.08556 ±0.0044 


4.78105 ±0.094 


xlO-i 




-1.21966 ±0.023 


-1.21143 ±0.0094 


-1.24077 ±0.035 


xlO"*^ 


% 


4.80077 ±0.51 


4.69457 ±0.0043 


4.30172 ± 1.1 


xlO-i 


% 


1.34418 ±0.081 


1.20476 ±0.064 


1.14697 ±0.21 


xlO"^ 


No cuts No ISR 



Table 2: Contributions to the total cross section for e^e hhW^W of the 
eight (non-)resonant channels of Ref. pO|, as obtained from the three algorithms 
documented in the text. For VEGAS, the setup ITMX = 5 and NCALL = 100000 
was adopted. 



The algorithms chosen for this test are Metropolis, VEGAS and RAMBO. In running 
RAMBO, we have not set up any special arrangement in dealing with the complicate 
phase space structure of the various channels |]T^^: we have let the algorithm 
generate four-momenta and weights and computed the estimate and error of the 
total cross section as described in the previous Subsection. As for VEGAS, we have 
adopted here the same mapping of the integration variables described in Ref, 



Concerning the Metropolis implementation, the expression in eq. (0) of Subsect. |2]l 



was used for computing the cross section. This method is applicable when the ME 
can be expected to have negative values. Also, the massive phase space volume was 
calculated beforehand and with insignificant errors. Notice that the A^mc statistics 
used here should be taken as representative of a value for which all three algorithms 
converge to the correct integrals (see Ref. []TU|,pU|). Indeed, if this is augmented. 
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errors diminish considerably in each case, though the relative performances among 
Metropolis, VEGAS and RAMBO remain basically unaffected. 



The eight terms %, i = 1, ...8, of eqs. (9)-(16) in Ref. pO|], integrated over the full 
four-particle phase space, are presented in Tab. 0, as obtained by using Metropolis, 
VEGAS (ITMX = 5 and NCALL = 100000) and RAMBO, with A^'mc = 500000. Here, 
it definitely is VEGAS to come out best, with second choice Metropolis and last 
RAMBO. The flaws of the latter should have been expected, as the algorithm is not 
adaptive so that it suffers from the presence of peaks rising over the phase space. 
This is particular evident in the Ty channel, which accounts for the very narrow 
H ^ bb resonance ||20| (recall that the Higgs width is just a few MeV)[3. In fact, 
when the resonant particles involved have a width of a few GeV (i.e., t, t, Z and 
W), such as in % with i = 2, ...6, the accuracy improves, unless two resonances 
have to be evaluated at once, those from top quark pairs in 7^. Here, the error 
does become very significant (about 10%). 

Metropolis behaves better than RAMBO, as its error is always smaller. It deals 
with single resonances rather satisfactorily for the quite low statistics used (when 
i = 4,7). No particular problems arise in Metropolis in dealing with negative 
weights (and rapidly changing interferences) either, i.e., % when i = 2,3,5,6,8, 
as here the typical error does not worsen in comparison to the cases in which the 
integrand function is definite positive (i.e., i = 4,7). However, Metropolis is no 
matching to VEGAS, particularly when multiple resonances are present, as in 7^. 

In the end, the careful mapping performed in VEGAS of all resonances and inter- 
ferences has paid off. However, one should notice the minimal involvement of 
the Metropolis implementation in this case, compared to the VEGAS one. in the 
Metropolis algorithm, there are no phase space Jacobian factors to be accounted 
for. On the other hand, we have stressed how they can efficiently be used in VEGAS 
to remove poles arising from the ME. Indeed, we will show in the next Subsection 
that, if an effort similar to that devoted to VEGAS here is employed for Metropo- 
lis as well (in 'teaching' to the algorithm the singular structure of the integrand 
function), then the performance of the latter can match that of the former. 



3.3 Variable energy, high multipHcity final states and cuts 

We dedicate this final Subsection of our numerical analysis to study the process 
gg —>■ bbti —>■ bbbbW^W~ bbbb]]l'^vi (where j represents any light-quark jet and 
jvt a lepton/neutrino), which was considered in Ref. as a QCD background 
to a possible charged Higgs discovery channel for the Large Hadron Collider (LHC). 
Numerical parameters and other inputs used for the runs were declared there. 



-^^Further notice that, being Mr < 2Mw, the Higgs decay channel H W^W in 7^ is not 
open. 
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In this test, there are three specific technical problems associated with the cal- 
culation of the total cross section. Firstly, the fact that the CM energy at par- 
tonic level is no longer a constant (contrary to the two previous examples): being 
a hadron-hadron process, two more integrations (in additions to those over the 
phase space) have to be performed, over the gluon momentum fractions, which 
evolve according to the PDFs inside the proton. Secondly, the very large number 
of final state particles, eight in total, which imposes a 21-dimensional integration 
over the phase spaceQ (beside the presence of various mass singularities, both in- 
frared and resonant poles). Thirdly, the introduction of severe reductions of the 
original integration region, as we have enforced in our simulations the same accep- 



tance and selection cuts recommended in Ref. [|2T|, the latter being implemented 
simply through theta functions^. 

The algorithms used in this example were only two: Metropolis and VEGAS. The lat- 
ter uses as usual a careful mapping around the heavy particle resonances, through 
the variable 

(j) = arctan I I (13) 

(where M and F represent the natural mass and width of the unstable particle, 
t,t or W^, with virtuality Q^), whose derivative is proportional to the resonant 
propagator itself: 

.<iQ^ (14) 



(g2 _ M2)2 + M2F2 



The setup of the former is as described in Subsect. p.3.3| . This is the most in- 
volved implementation of the Metropolis algorithm which was tested. In this case, 
one needs an approximate ME with known cross section at a certain energy: this 
was constructed by simply using the two W- and the two t-resonances. There- 
fore, the Metropolis approach adopted here can be viewed as the equivalent of the 
VEGAS mapping enforced through eqs. ([T^)-(|T^. The cross section associated to 
this auxiliary ME at one energy was calculated numerically beforehand. In the 
following, we assume the latter to be known with arbitrary small error. In this re- 
spect, we should also mention that in our actual ME we have ignored interference 
effects between the two subsets of Feynman diagrams that only differ in the ex- 
change of the four-momenta and spins between the two 6-quarks (or, equivalently, 
the two 6-antiquarks) in the final state, because of their indistinguishibility in the 
detector, and a minus sign, because of the Fermi-Dirac statistics (in other terms, 
the Pauli principle). In fact, their effects on the total and differential cross sec- 
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One degree of freedom being absorbed int o the flat azimuthal integration about the incoming 



beam direction (see discussion in Subsect. -iA), and already accounting for the PDF convolution. 

^^In fact, it turns out impossible to map the entire phase space in terms of the kinematic 
quantities whose range is being cut, and not any more efficient to use only one or two of these as 
integration variables. 
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tions are negligible [?T|. Besides, their integration would pose further, unnecessary 
complications. 



^a/a{gg ^X^ bbbbjji^u,) (%) 


y/§ (TeV) 


Metropolis 


VEGAS 


X only top-antitop radiation, no PDFs 


0.6 


1.92 


0.48 


1.0 


4.25 


0.99 


1.4 


6.88 


1.52 


X full ME, no PDFs 


0.6 


0.77 


0.32 


1.0 


2.15 


0.61 


1.4 


2.96 


0.72 


X full ME, PDFs 


Vrs 


3.58 


1.53 


Kinematics cuts from Ref. 


21] 



Table 3: The relative error on several cross sections associated to the process 
gg ^ X ^ bbbbiii^Ui at y/s = 14 TeV, as obtained from the two algorithms 
documented in the text, each using about 10^ MC points (all passing the default 
cuts of Ref. |]2T|). We have verified that actual cross sections (not shown here) are 
statistically consistent between the two algorithms. 

As we have already digressed to some length about the relative ability of the two 
algorithms to adapt, we make our primary concern in this test that of comparing 
the size of the errors associated to the integrals in each case. In order to render 
the comparison consistent, regardless of the actual value of point generated, Nmc, 
we always compute the integrals for a given statistics in both cases, e.g., 10^ (that 
is, the latter is the approximate number of MC points that actually pass the cuts). 
We proceed to obtaining the final result by steps, in order to assess whether one 
algorithm outperforms the other in some specific task. We start by isolating a 
gauge-invariant substructure of the original ME, only comprising those diagrams 
(eight in total) in which the off-shell gluon, g*, eventually splitting into bb pairs, is 
emitted by either of the top (anti) quarks, the latter finally decaying semileptoni- 
cally. Moreover, we fix the partonic CM energy, i.e, = constant, thus removing 
for the time being the integrations over the gluon PDFs. As the corresponding 
cross section has little meaning physics-wise, we only plot the relative error as 
obtained from the two algorithms, e.g., at y/l = 600, 1000 and 1400 GeV. The 
upper part of Tab. |^ shows the rates for the simpler subprocess just described, i.e., 
gg ^tt ^ g*tt bbtt bbbbW^W- bbbbHi^ue. At all energies, VEGAS yields 
a smaller error than Metropolis, by about a factor of four. Although Metropo- 
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lis is still outperformed by VEGAS in the size of the relative error of the various 
integrations, one should notice that the differences between the two algorithms 
have diminished substantially: compare to the rates in Tab. 0. In both cases, the 
relative error increases with In order to understand this effect, it should be 
recalled that, although the final state particles are jets and leptons (thus with neg- 
ligible rest mass as compared to the values used), the two (anti)top resonances 
involved impose that the cross section would drastically fall to negligible levels if 
~ 2(mt + mb). Effectively, the volume of the phase space associated with the in- 
tegration performed at y/S = 600 GeV is much smaller that that spanned when y/I 
is 1400 GeV, where the 8-particle phase space can stretch much further away from 
the bbti threshold at 360 GeV. Therefore, one would conclude that the tendency 
of both integrators to giving smaller errors at lower CM energies is a consequence 
of the fact that the integrand function fluctuates more at larger 



As a second step of our test, we have introduced the full ME, see Ref. in place 
of the reduced one considered so far, where by 'full' we intend the one obtained 
by allowing for the attachments of the off-shell g* — > bb current to the gluon lines 
too, but with the process still proceeding via the production of two top (anti)top 
quarks, and without the mentioned interferences. This 2^8 ME consists of 36 
basic Feynman diagrams. Not to complicate things further, we again express the 
gluon PDFs through 5-functions centered around one and fix the partonic energy 
y/I at the usual three values: 600, 1000 and 1400 GeV. Results are presented in 
the middle part of Tab. ^ The relative performances of the two algorithms are 
rather similar to the previous case, though the overall error has diminished at each 
energy in both Metropolis and VEGAS. This effect is presumably the consequence of 
the fact that the additional contributions to the ME fill regions of the phase space 
previously empty, these smoothly interpolating into those typical of the kinematics 
of gluon emission from (anti)top quarks. 

The final step is the integration over the full ME, including the convolution with 
the PDFs, i.e., s = X1X2S = ts, with y/s = 14 TeV. In VEGAS, the latter was 
done by adding the two integrations over the two gluon momentum fractions, i.e., 
Xi and X2, (so that NDIM = 21) and the call to the PDF numerical package. 
In Metropolis, we have proceeded as follows. One of the integrations over the 
momentum fractions was performed beforehand. This was done by first changing 
the momentum fraction variables into the logarithm of the (squared) CM energy 
and rapidity of the two gluons (as in VEGAS): 

dxidx2 , , , ... /x7 

= aln[XiX2)dln . — . 

Xi X2 V ^2 

The rapidity spectrum was then integrated at fixed CM energy, this yielding a 
distribution in the latter variable only. This was in turn convoluted into the weight- 
function. The complete integration was finally done by combining the variable 
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energy and parallel integration methods, as described in Subsect. p^.3.3| . The region 
400 GeV ~ ~ 2 TeV was chosen, large enough to gather the main contributions 
to the cross section. From the lower part of Tab. ^, one can see that, even in 
presence of the complete differential structure of the cross section, VEGAS performs 
better than Metropolis in terms of accuracy, but the relative size of the error has 
gone down, to a factor slightly larger than two only. The absolute size is larger 
than in the previous test for both algorithms, a consequence of the additional 
dependence upon the gluon PDFs. 

As intimated at the end of the previous Subsection, if a similar sort of care is 
devoted to both algorithms in order to better account for the singular structure of 
the MEs, one should conclude that, although VEGAS is still better in minimising 
the accuracy of the integration. Metropolis represents at least a viable alternative. 
This becomes particularly true if one finally considers that, even in this context. 
Metropolis displayed a better tendency than VEGAS to adapt in high dimensions 
and over a reduced phase space, as already seen elsewhere. In other words, the 
same accuracy can be achieved by the former with much less CPU-time effort than 
for the latter. 



4 Conclusions and outlook 

Adaptive MC programs have become an indispensable tool in high energy parti- 
cle physics, in order to calculate reliably decay and scattering cross sections over 
multi-dimensional phase spaces. Indeed, as the energy reach of modern particle 
accelerators grows larger, both the number of particles that can be accommodated 
in the final state and that of the channels through which they can be produced, 
increase rapidly. Under these circumstances, it is evident that simple generalisa- 
tions of well known one-dimensional integration methods are no longer applicable, 
given the huge number of points that would be needed in order to overcome the 
high complexity of the integrand functions (let alone the use of analytic methods, 
if one further considers the need of accounting for non-trivial reductions of the 
integration volume, because of unavoidable experimental cuts). 

Although several algorithms already exist on the market nowadays, which perform 
their task sufficiently well to have enabled stringent tests between theory and 
experiments, it is clear that the availability of new implementations is a need never 
exhausted: if only for cross-checking purposes. For this and other reasons, we have 
developed several new implementations of an old MC technique: the so-called 
Metropolis method. Although being another adaptive MC algorithm, it can boast 
at least one radically different feature with respect to traditional approaches: in the 
latter, the integration over the phase space takes place over a grid; in the former, it 
evolves along a path. This discretisation of the integration volume has historically 
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appealed to the solution of problems in statistical mechanics and Lattice Gauge 
Theories. However, there is basically no reason why such a technique can not 
be applied successfully also to the four-momentum phase space of Quantum Field 
Theories. 

As a matter of fact, we have here demonstrated the high potential of the Metropo- 
lis method in dealing with realistic problems arising in modern particle physics 
phenomenology. It not only fulfills the basic criteria of accuracy and efficiency 
required to any algorithm, no matter the number of dimensions involved in the 
integration, but it has also been shown to outperform in some cases other, already 
widely diffused MC programs. Moreover, being its speed matter of no concern 
at all, the Metropolis algorithm could well be used in (parton level) MC event 
generators too. 

There can certainly be drawbacks in the application of the method. The most im- 
portant being the need at times of optimising the implementation of the Metropolis 
algorithm to solve a specific problem, which can require prior knowledge of the be- 
haviour of the integrand function. But even then, the reader should acknowledge 
that it is becoming more and more rare in numerical simulations that one can 
afford to rely solely on the ability of whatever algorithm in adapting itself to such 
subtle effects as interferences, finite particle widths, irreducible backgrounds, etc., 
as some of our examples should clearly have demonstrated. 

We have encoded the various implementations of the Metropolis algorithm dis- 
cussed in this paper in a C++ program, that we make available to the public upon 
request. This code makes extensive use of the CLHEP classes |]2^ for handling 
Lorentz four-momenta. 

Before closing, we will describe two possible improvements that could be considered 
in the future in order to further increase the efficiency of the Metropolis algorithm 
in QFT. Both these methods are widely used for conventional statistical physics 
problems, and shown to be necessary tools in many cases. 

The first method relies on the introduction of a fictitious reciprocal temperature, j3. 
This method has seen two different implementations called 'simulated tempering' 
| n| and 'parallel tempering' [^0. In our case, the phase space would be enlarged 



with (3 G [0, 1] as a new parameter, inducing a modified weight-function, 

w{x^ (3) = w^{x). 

Clearly, the integration is simpler with (3 close to zero, while we want it instead 
to be performed at /? = 1. This way of enlarging the phase space allows for round 
tours into the low-/5 region where the configuration can be updated more freely. 



^^The difference between the two lies essentially in the way the configuration space is enlarged 
by means of (3. 
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and thus the correlation length is shorter. This is very much in the same manner 
as when the total energy, E^ot, was used as a free parameter (see Subsect. |2.3.2| of 



this article). The advantage of instead using (3 is that the efficiency of the method 
would not be as much dependent on the energy variations of the cross section. In 
any case, this method should, and could rather easily, be tested also in the case of 
integrations over a four-momentum phase space. 

The other technique that we want to mention is known as the 'Hybrid Monte Carlo' 
method In this case, derivatives of the weight function are used to move large 
distances in configuration space, but still avoiding the un-favoured regions where 
the weight function is small. In many cases, this method has increased the efficiency 
of the Metropolis algorithm drastically and could very well be unavoidable if much 
more numerically demanding problems are to be treated in the QFTs of the future. 
Unfortunately, derivatives of more complicated matrix elements are not available 
today. It is also at present not clear to the authors how this hybrid MC should be 
implemented to perform integrations over a four-momentum phase space. 

A final, more general remark, is the following. Throughout this article, we have 
described several implementations of the Metropolis algorithm, all using a stepping 
procedure designed in such a way that the final state particles are kept on mass- 
shell. An alternative approach would be to map their four-momenta onto a set 
of independent variables and perform the random walk in the new phase space. 
This way, one could benefit from the use of mappings which can cancel the ME 
singularities (e.g., see eqs. (|^) and (p^if ) in Sect. |3.3|) . A possible implementation 
of this last strategy is now under consideration. 



Appendix: a proof for restricted phase spaces 

In this Appendix, we will show that, in case the integration is to be done in a 
subspace of the original phase space, the PS^-ed of Subsect. |2.1|^°|, the Metropolis 



condition (a) is fulfilled, if for any two points inside PS'red, there exists a continuous 
curve r, which connects the two points and is contained in PSj-cd- 

For every point x along F, there is an open disc D{x, r), with radiu^ r, which is 
wholly contained in PS^ed- First, we show that we can choose a. 6, < 6 <r, so 
that all points on the surface {y; \x — y\ = 6} can be reached from x in n — 1 steps 
which are all inside D. 

Suppose, as before, that x = {pi} and y = {p-} and let Aj = p'i — Pi- We can now 
move from xtoyinn — l steps, by transferring the Aj one by one to Pn for z = 1 to 
i = n — 1. This will take us through n — 1 points: a; = a;i = y. 

'^^Which we assume here to be an open subspace. 

^^With the metric of the EucUdian space of the three-momenta. 
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Let d be the largest of |x — Xj| in this path. Choose an a such that < a < r /d 
and get a new 6' = a6. Then the point y' = x + a{y — x) lies on the corresponding 
surface and can be reached by a path inside the disc D. If y is the point on 
{y; \x — y\ = 6} which gives the largest d = d and with 0<d<r/d, we have that all 
points in S* = {y; \x — y\ = d6} can be reached in a finite number of steps, which 
are all contained in D. 

A point where F crosses 5* has a finite distance {d6) to x and it can then be reached 
with a finite number of steps, all inside PSj-ed- In this way, starting from any point 
in PSj-ed, any other point can be reached with a finite number of steps inside PS^ed, 
by taking finite steps along the corresponding F. This concludes our proof. 
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Abstract 



A proof presented in the original paper was incorrect. We outline here an 
alternative procedure. 



The new proof 



One of the conditions required to the stepping procedure implemented in the 
Metropohs algorithm, as described in Ref. (point (a) in Sect. 2.1), was that 
any point {p[} in the n-particle phase space can be reached from any other point 
{pi} in a Gnite number of steps. The proof in the main article relies on the false 
statement that, for any two on-mass-shell particles with four-momenta pi and p2, 
the whole of the surface of conserved pi + p2 can be reached by the identical par- 
ticles by rotating them in their centre-of-mass (CM) frame. This is however not 
true, because it would in general be possible to transfer a four-momentum which 
conserves pi +p2 but puts one or both particles out of the mass-shell. (The rotation 
procedure itself, in contrast, does preserve the on-mass-shell condition). Below we 
present a correct proof of the above statement. 

It is clear that the whole of the two particle phase space can be reached in a 
single step, so let us assume that the number of particles (n) is at least three. We 
will later on show that it is possible to move any particle (e.g., the n:th.), from 
its original four-momentum pn to the final one in a finite number of steps, by 
exchanging some momentum with the other particles. This procedure will change 
the other four- momenta to new values {Pi}- The rest of the problem is then to 
move these to their final values {p-}. For (n — 1) = 2, this can be accomplished in 
one step. For n > 3, we can move the (n — l):th four-momentum to its final value 
p'n-i, and so on. 

Let us first discuss what four-momenta are possible for a given particle, in the case 
of unrestricted phase space. Consider again the n:th particle, and let s be the 
total invariant energy squared, be the energy squared of the system of all the 
other particles and let m„ be the mass of the n:th particle. In the CM frame, the 
magnitude of the three-momentum of the n:th particle is given byQ 




with 

A(a, b, c) = a? + 1)^ + c^ - 2ab - 2bc - 2ca, 

and it can have any direction. Depending on the value of s„_i, P„ can have any 
value between (for s„_i = (-\/s — m„)^) and its maximum, P™^^, which occurs at 
s„_i = (J27=i ''^i)'^ (^h^ minimal value of s„_i). 

The four-momentum p„ is, from its mass-shell condition, restricted to a three- 
dimensional sub-space. We will now show that, in case P„ is smaller that its 
maximal value, there exists a three-dimensional on-mass-shell neighbourhood of 

^^Hereafter, we use P to indicate the magnitude of a momentum in the CM frame and p to 
denote a four-vector in any frame. 
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Pn in which any point can be reached in two steps, by exchanging momentum with 
two other particles (denoted below by pi and ^2)- 

Let ki denote the exchanged four-momentum from pi to Pn- Our condition for 
the exchanged momentum is that it keeps both of the particles on their respective 
mass-shell. Consequently, with Pn = Pn + ki we have = = + A;^ + 2kiPn 
kf — —2kiPn- Similarly, we have the condition kf — 2kipi, in order to keep pi on 
the mass-shell. These two conditions restrict ki to lie on a two-dimensional surface 
(corresponding to the two rotational degrees of freedom in the CM frame of Pn 
and pi.) We will now show that two particles pi and p2 can be chosen so that the 
corresponding infinitesimal four-momentum transfers eki and ek2 (e 0) to Pn are 
such that k = eki + ek2 has three degrees of freedom. This implies that Pn + k (for 
a finite k) covers an open on-mass-shell neighbourhood of p„. In the limit e — > 0, 
the on-shell conditions impose 



Under such conditions, the four-vectors ki are restricted to lie in two-dimensional 
planes. The vector k will then have more than two degrees of freedom, unless the 
planes for ki and k2 coincide. This would in turn imply that the four-momenta 
Pi and P2 are linearly dependent, that is, one can write pi = ap2 for a number a. 
Consequently, it is enough to show that there exist two other particles such that 
their four-momenta are not proportional. This condition is equivalent to say that 
Pn, as was assumed before, is smaller than its maximal value. 

To demonstrate this statement, we will show that reaches its minimum if and 
only if all n — 1 four-momenta are proportional to each other, that is, there are 
numbers aij so that one can write pi = aijpj, for any i,j < n — 1. This is well 
known in case all particles are mass-less (s„_i = 0). In case any of them has a 
mass, one can boost the entire system to the rest frame of such a particle. The 
minimal value of s„_i is then reached if all particles have zero three- momentum in 
the new frame. 

This completes the proof that, in case P„ < P™^, there is a three-dimensional 
neighbourhood of p„ which can be reached in two steps. From this, wc draw the 
conclusion that the final point p'^ can be reached in a finite number of steps, in case 
both Pn and P'^ arc smaller than the maximum. If this was not the case, then there 
would exist a four-momentum p" , with P" < P^^, which p^ can come arbitrarily 
close to, but never reach. But this contradicts the fact that a neighbourhood of p" 
is reachable from p" (and vice versa). 

The final step is to show that it is always possible (in case n > 2) to move from 



Pnk2 

Piki 

P2k2 



0, 
0, 
0, 
0. 
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a point with P„ = P^^"^ to another point with < P^^^. This follows from the 
condition pi = aijpj, that leaves one degree of freedom for the four-momentum of 
one of the particles, for fixed values of the others. Since the four-momentum trans- 
fer k has two degrees of freedom, one can, in one step, violate the proportionality 
condition above, with the consequence that P^ < P^^. 
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